function [bestChromosomeCell, bestInternalFitness, fitnessHistoryInternal] = ImprovedGeneticAlgorithm(DSM, params)
%IMPROVEDGENETICALGORITHM 实现改进的 GA 以实现模块化。
%
%   Inputs:
%       DSM    - M x M 设计结构矩阵
%       params - 包含算法参数的结构
%
%   Outputs:
%       bestChromosomeCell     - 包含最佳解决方案矩阵的单元格数组。
%       bestInternalFitness    - GA 的内部最佳适应度（最大化，即 -min_objective）。
%       fitnessHistoryInternal - GA 每代内部最佳适应度值的向量。

    M = params.M;

    fprintf('Initializing population...\n');
    Population = initializePopulation(M, params,DSM); 
    fprintf('Population initialized with %d individuals.\n', length(Population));

    bestChromosomeCell = {}; % 存储为潜在可变大小矩阵的单元格
    bestInternalFitness = -inf; % GA 最大化了这一点
    fitnessHistoryInternal = zeros(params.MaxGenerations, 1);

    for gen = 1:params.MaxGenerations
        fprintf('\n--- Generation %d/%d ---\n', gen, params.MaxGenerations);

        fitnessValues = zeros(length(Population), 1);
        for i = 1:length(Population)
            if isempty(Population{i}) % 安全检查
                fitnessValues(i) = -1e8; % 空染色体的适应度非常差
                continue;
            end
            % 计算适应度值 (calculateFitness) 返回实际目标的负值
            fitnessValues(i) = calculateFitness(Population{i}, DSM);
        end

        [currentGenBestInternalFitness, currentGenBestIdx] = max(fitnessValues); % 数组 fitnessValues 中找出当前代
        % （generation）中最优的适应度值以及其对应的索引位置
        if currentGenBestInternalFitness > bestInternalFitness
            bestInternalFitness = currentGenBestInternalFitness;
            bestChromosomeCell = {Population{currentGenBestIdx}}; % 将矩阵存储在单元格中
            fprintf('New GA best internal fitness: %.4f (True Min Objective: %.4f)\n', ...
                    bestInternalFitness, -bestInternalFitness); % Display both
        end
        fitnessHistoryInternal(gen) = bestInternalFitness;

        [~, sortedIdx] = sort(fitnessValues, 'descend'); % 将 fitnessValues 按从大到小排序，并获取其排序后的索引顺序
        elites = cell(params.EliteSize, 1);
        eliteIndices = sortedIdx(1:params.EliteSize);
        for i = 1:params.EliteSize
            elites{i} = Population{eliteIndices(i)};
        end

        fprintf('Clustering population...\n');
        subPopIndices = hierarchicalClustering(Population, fitnessValues, params);
        numSubPops = length(subPopIndices);
        fprintf('Clustered into %d sub-populations.\n', numSubPops);

        newPopulation = cell(params.PopulationSize, 1);
        newPopCount = 0;

        for i = 1:params.EliteSize
             newPopCount = newPopCount + 1;
             newPopulation{newPopCount} = elites{i};
        end

        offspringNeeded = params.PopulationSize - params.EliteSize;
        offspringGenerated = 0;
        avgFitness = mean(fitnessValues(isfinite(fitnessValues))); % Exclude Inf/-Inf for avg

        subPopAvgFitness = zeros(numSubPops, 1);
        for sp = 1:numSubPops
           valid_indices_in_subpop = subPopIndices{sp}(isfinite(fitnessValues(subPopIndices{sp})));
           if ~isempty(valid_indices_in_subpop)
               subPopAvgFitness(sp) = mean(fitnessValues(valid_indices_in_subpop));
           else
               subPopAvgFitness(sp) = -inf; % Or some other default for empty/all-Inf subpops
           end
        end
        
        offspringPool = cell(offspringNeeded * 2, 1); % Allocate more initially for pairs
        tempOffspringCount = 0;

        while tempOffspringCount < offspringNeeded * 2 % Generate pairs for crossover
            if offspringNeeded == 0; break; end % No offspring if pop is all elites

            parentIdx1 = tournamentSelection(fitnessValues, params.TournamentSize);
            parentIdx2 = tournamentSelection(fitnessValues, params.TournamentSize);
            Parent1 = Population{parentIdx1};
            Parent2 = Population{parentIdx2};

            Child1 = Parent1; 
            Child2 = Parent2;

            if rand < params.CrossoverRate
                [Child1, Child2] = adaptiveCrossover(Parent1, Parent2, params);
            end

            Child1 = adaptiveMutation(Child1, DSM, fitnessValues, avgFitness, params);
            Child2 = adaptiveMutation(Child2, DSM, fitnessValues, avgFitness, params);
            
            Child1 = variableNeighborhoodSearch(Child1, DSM, params);
            Child2 = variableNeighborhoodSearch(Child2, DSM, params);
            
            if ~isempty(Child1)
                tempOffspringCount = tempOffspringCount + 1;
                offspringPool{tempOffspringCount} = Child1;
            end
            if tempOffspringCount < offspringNeeded * 2 && ~isempty(Child2)
                tempOffspringCount = tempOffspringCount + 1;
                offspringPool{tempOffspringCount} = Child2;
            end
            if tempOffspringCount >= offspringNeeded * 2; break; end
        end

        % Fill the new population with generated offspring
        for i = 1:min(tempOffspringCount, offspringNeeded) % Take only needed amount
             newPopCount = newPopCount + 1;
             if newPopCount > params.PopulationSize; break; end 
             newPopulation{newPopCount} = offspringPool{i};
        end
        
        % If not enough offspring were generated (e.g. due to errors), fill with random
        while newPopCount < params.PopulationSize
            newPopCount = newPopCount + 1;
            num_rand_modules = max(2, round(sqrt(M)));
            newPopulation{newPopCount} = generateSingleChromosome(M, num_rand_modules, 'random');
        end
        
        Population = newPopulation(1:params.PopulationSize); % Ensure exact size

        if mod(gen, params.AprioriFrequencyTheta) == 0 && gen > 0 && ~isempty(Population)
            fprintf('Running periodic Apriori and Improvement...\n');
            try 
                currentFitnessValuesApriori = zeros(length(Population), 1);
                 for i_ap = 1:length(Population)
                     if ~isempty(Population{i_ap})
                        currentFitnessValuesApriori(i_ap) = calculateFitness(Population{i_ap}, DSM);
                     else
                        currentFitnessValuesApriori(i_ap) = -1e8; % Penalty
                     end
                 end
                [~, sortedIdxMining] = sort(currentFitnessValuesApriori, 'descend');
                numEliteMining = max(1, min(length(sortedIdxMining), floor(0.1 * params.PopulationSize))); 
                eliteIndicesMining = sortedIdxMining(1:numEliteMining);
                
                if ~isempty(eliteIndicesMining)
                    transactions = convertToTransactions(Population(eliteIndicesMining));
                    frequentPatterns = aprioriMiner(transactions, params.AprioriMinSupport); 

                    if ~isempty(frequentPatterns)
                        fprintf('Found %d frequent patterns. Improving...\n', length(frequentPatterns));
                        [~, sortedCurrentIdxWorst] = sort(currentFitnessValuesApriori, 'ascend'); % Find worst
                        numToImprove = max(1,min(length(sortedCurrentIdxWorst), floor(params.PopulationSize * 0.1)));
                        indicesToImprove = sortedCurrentIdxWorst(1:numToImprove);
                        Population = improveWorstSubpopulation(Population, indicesToImprove, frequentPatterns, params);
                    else
                        fprintf('No frequent patterns found.\n');
                    end
                end
            catch ME_apriori
                 warning('Error during Apriori/Improvement: %s', ME_apriori.message);
                 if strcmp(ME_apriori.identifier, 'MATLAB:badsubscript')
                     fprintf('Bad subscript error likely due to empty population or indexing issue.\n');
                 end
            end
        end
    end % End main generation loop

    fprintf('\nGA Finished.\n');
    fprintf('Best GA Internal Fitness (Maximized): %.4f\n', bestInternalFitness);
    fprintf('Corresponding True Best Minimum Objective: %.4f\n', -bestInternalFitness);

end % End main function ImprovedGeneticAlgorithm